function A = sparse_matrix_function(injndiag,ijndiag,ipjndiag,injdiag,ijdiag,...
        ipjdiag,injpdiag,ijpdiag,ipjpdiag,I,J)
    
    % V_{0,j} = V_{1,j}    
    ind = (1:I:size(ijndiag,1));
    ijndiag(ind) = ijndiag(ind)+injndiag(ind);
    ijdiag(ind) = ijdiag(ind)+injdiag(ind);
    ijpdiag(ind) = ijpdiag(ind)+injpdiag(ind);
    injndiag(ind) = 0; injdiag(ind) = 0; injpdiag(ind) = 0;
    
    % V_{I+1,j} = V{I,j}
    ind = (I:I:size(ijndiag,1));
    ijndiag(ind) = ijndiag(ind)+ipjndiag(ind);
    ijdiag(ind) = ijdiag(ind)+ipjdiag(ind);
    ijpdiag(ind) = ijpdiag(ind)+ipjpdiag(ind);    
    ipjndiag(ind) = 0; ipjdiag(ind) = 0; ipjpdiag(ind) = 0;
    
    % Adjust along jth dim now
    ind = 1:1:I; % V_{i,0} = V_{i,1}
    ipjdiag(ind) = ipjdiag(ind)+ipjndiag(ind); % upper diag of center
    ijdiag(ind) = ijdiag(ind)+ijndiag(ind); % main diag of center
    injdiag(ind) = injdiag(ind)+injndiag(ind); % lower diag of center
    ipjndiag(ind)=0; ijndiag(ind)=0; injndiag(ind)=0;

    ind = (I*J-I+1):1:I*J; % V_{i,J+1} = V_{i,J}
    ipjdiag(ind) = ipjdiag(ind)+ipjpdiag(ind);
    ijdiag(ind) = ijdiag(ind)+ijpdiag(ind);
    injdiag(ind) = injdiag(ind)+injpdiag(ind);
    ipjpdiag(ind)=0; ijpdiag(ind)=0; injpdiag(ind)=0;

    % Adjust coefficients to fit sparse matrix (particularity of spdiags)
    ijndiag = [ijndiag(I+1:end); zeros(I,1)];
    injndiag = [injndiag(I+2:end); zeros(I+1,1)];    
    ipjndiag = [ipjndiag(I:end); zeros(I-1,1)];
    % V_{i,j} coefficient: use full diagonal
    injdiag = [injdiag(2:end); 0];
    ipjdiag = [0; ipjdiag(1:end-1)];
    ijpdiag = [zeros(I,1); ijpdiag(1:end-I)];
    injpdiag = [zeros(I-1,1); injpdiag(1:end-I+1); 0]; 
    ipjpdiag = [zeros(I+1,1); ipjpdiag(1:end-I-1)];
    
    % Sparse representation
    A = spdiags(ipjpdiag,I+1,I*J,I*J) + ...
        spdiags(ijpdiag,I,I*J,I*J) + ...
        spdiags(injpdiag,I-1,I*J,I*J) + ...
        spdiags(ipjdiag,1,I*J,I*J) + ...
        spdiags(ijdiag,0,I*J,I*J) + ...
        spdiags(injdiag,-1,I*J,I*J) + ...
        spdiags(ipjndiag,-I+1,I*J,I*J) + ...
        spdiags(ijndiag,-I,I*J,I*J) + ...
        spdiags(injndiag,-I-1,I*J,I*J);    
    